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Abstract: — A quantitative .prediction of the thermal field and. the location of the melt-crystal interface .. 
.: during growth requires a precise knowledge of the heat transfer taking place ,in the entire furnace. The^ 
problem is highly complex since it* involves an accurate' calculation 1 of radiation. between the different 
surfaces' and conduction in the constituents. Radiative exchanges are calculated with the assumption of : 
dirruse^surfaces : ^and with use of." a viewed and hidden part . algorithm : together with a Galerkin 
discretization. The model has been- extended for taking semi-transparent materials into account The shape . • 
of the liquid-^sqlid mterface is a variable of the problem and is calculated as being the melting isotherm: 
Examples of germanium and gallium arsenide furnaces are analysed, showing the efficiency -of the method." ' 



1. INTRODUCTION 

Most semiconductor crystals are grown in Czo- 
chralski pullers, in 'which a cylindrical crystal is pulled 
from the suiface ofithe melt. Gallium arsenide and 
indium phosphide crystals axe usually grown with the 
LEC (liquid-encapsulated Czochralski) process. The 
calculation of heat transferrin "such furnaces is complex 
in view of the strong coupling between diffusive, 
advective and radiative exchanges. A review of the 
various aspects affecting heat transfer in a Czochralski 
furnace may be found in ref. [1]. Over the last four 
years; a major effort has been accomplished for con- 
sidering the whole furnace as a system where the 
growth of the crystal is simulated on the basis of a 
reduced number of controllable parameters, such as 
the power input in the heater^ the coolant temperature,' 
the pulling rate and the diameter of the crystal. In a 
series of papers, Derby etidl.*[2-6\ developed a sim- 
plified model of the furnace radiating towards ah 
ambient temperature. They studied the -dynamics of 
the growth process and its control, and considered the 
variation of the crystal radius during growth. A more 
elaborate thermal^ into account radiative 

exchanges between the various constituents was later 
elaborated by % 'Athertpn : e? al.. [7], but the ambient 
temperature on tie" outer walls was still part of the 
data. In an earlier paper, Srivastava et al. [8] had also 
used a radiative model on a reduced geometry, with 
an analytical ^calculation of- the view factors. More 
recently, Mbtakef and Witt [9] and Motakef [10, 11] 
developed a thermal model which is very similar to 
that in ref. P), except that the radius of the crystal is 
uniform:, -vi:;".: 

In these various papers, a certain .degree, of inde- 
terminacy is left on the distant surfaces such as the 
wall of the furnace or' the bottom plate; the effect of 



including 1 additional components in the heater, such 
as reflectors above trie crucible,, is not straightforward. 

In the present paper, we wish to pursue .the analysis 
initiated by Wputers. [12]. The k method has briefly 
been reviewed in ref.. [13], where we described a' heat 
transfer model of the furnace which is self ^contained ; 
mput parameters other than the. coolant temperature 
and the power input — or. the ..pulling rate— are not 
required, provided the complete geometry (including 
the. uniform radius .of the crystal) is . fully charac- 
terized. Since the -publication of ref s. .[12, 13], the 
method has been refined and several practical appli- 
cations have been published [14-17]. In. what follows, 
we wish to. give a.complete description of the thermal 
v model and the numerical, method used in these papers. 

A major ingredient of the /model is the calculation 
of the.radiatiyeh transfer... ^^Aj^yti.cal.metho,<k. were 
eliminated a priori in view of the geometrical com- 
plexity: of common furnaces and. the nonlinearity of 
the problem. A first numerical approach would be the 
Monte Carlo method, which lies on a. statistical basis 
[18, 19]. Such : a method is well adapted for solving 
complex problems,; but.it has not been retained in view 
of its .high cost for a sufficient accuracy. A .second 
approach is the net-radiation method ; when the latter 
is. combined with zonal analysis .[20V23],. the surface 
of the enclosure is divided into a number of isothermal 
patches, while radiative .transfer is. calculated between 
every .pair of such patches. . First, one calculates the 
configuration (or view) factor associated with every 
pair, while taking into account the viewed and.hidden 
parts of the enclosure. Matrix operations then provide 
the., total-exchange areas , : (or .the Gebhart factors) 
between the pairs; this means that successive reflec- 
tions are now taken into account. The total-rexchange 
areas produce the relationships between net fluxes and 
radiative emitted -fluxes, i.e. the matrix relationship 
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between net fluxes and fourth powers of temperatures 
on the enclosure. Such an approach was followed by 
Atherton et al. [7], where circular cross-sections of 
the crucible are replaced by octagons. . _ 

In our work, we have taken full account of the 
circular symmetry of the furnace. The temperature 
field on the enclosure is represented by means of quad- 
ratic one-dimensional finite . elements. . Qalerkin's 
method [24] is usee! for solving the radiative integral ' 
equation ; it leads to a high accuracy which is closely 
related to that of the numerical integration. Let us 
in particular emphasize the importance of radiative 
exchange near the tri-junction, where the crystal meets . 
the melt View factors between infinitesimal conical 
surfaces are calculated analytically by of an 

original method which is equivalent to but differs from 
Garot and Gendre's method. [25]. A major advantage 
of our algorithm is that the full analysis.is performed 
in a planar cross-section of the furnace, and -is thus 
an appropriate tool for furnace design. ? , 

A study of the LEC growth requires to take into 
account the presence of a semi-transparent oxide 
layer. While the limit cases of transparency and opac- 
ity are^studied in ref. [10], we have applied the band- 
energynlethod [19], considering that the layer is trans- 
parent for some ranges of wavelengths and opaque 
for ' others. : We will find' that tie "transparency 
coefficient has important implications upon 5 the tem- 
perature distribution and the thermal stresses. 

•A significant problem of using semiconductors for 
high quality electronic substrates is the presence of 
dislocations. They . .are usually related to thermal 
stresses [26, : 27] generated during growth : Which, 
beyond a critical level, give rise to plastic -defor- 
mations. Thus, a common way of estimating the den- 
sity of dislocations, at least on a qualitative basis, 
has been to calculate thermal stresses during growth. 
Analytical methods have been used by Kobayashi and 
Iwaki [28] and Jordan' er al. [29, ;30] v While Duseaux 
[31] used a finite element method for this purpose. In 
such studies, the thermal boundary conditions oh the 
crystal are selected at the outsetand take little account 
of the surrounding ambience and of the shape of the 
interface.' An 'examination of thermal stress contours 
shows however that an accurate evaluation of such 
stresses requires in turn an accurate description of 
the temperature field. Methods based 'on measured 
temperatures [32] lose the predictive power which is 
desired for modifying the geometry of the furnace. An 
analysis of the thermal stresses on the basis of a global 
heat' transfer analysis has been developed by Motakef 
and Witt [9] and Motakef [10, 11]. Similarly, a finite 
element prediction of thermal stresses is coupled to 
our global heat transfer calculation. For this reason, 
our model is a powerful tool for acting on the mag- 
nitude of the temperature gradients in order to 
decrease the dislocation density (or the thermal 
stresses) .[14, 16]. V 7 

'At this' stage, we have not mentioned the motion of 
the meif caused by natural and forced convection. 



Although advective phenomena have been taken into 
account in our earlier publications [15, 17, 33], we 
have not included their description in the present 
paper for the sake of brevity., It should be recalled 
■that the low Pfandtl number characterizing metallic 
melts reduces the influence of convection on the tem- 
perature field, particularly in the case of germanium 
growth. 

* Since our goal here is to present the method for 
calculating heat transfer in the furnace, the paper 
is organized along the basic steps ' of the algorithm. 
"Having introduced definitions and basic hypotheses 
in Section 2, we present our algorithm for radiative 
enclosures in Section 3. Heat conduction is studied in 
Section 4, while the solid-liquid interface is e xamin ed 
in Section 5. In Section 6 we : elaborate; oritlie impor- 
tant topic of selecting a strategy for solving the global 
system of non-linear equations ; in particular, we dis- 
cuss whether one should impose the power input or 
foe pulhng rate. Having briefly recalled the i calculation 
of thermal stresses in Section 7, we analyse two exam- 
ples in Section 8, related to the growth of germanium 
and gallium arsenide crystals. 

2. BASIC HYPOTHESES AND DEFINITIONS 

The geometry of the global problem that' we wish 
to solve isillustrated in Fig: l, : where the components 
of a typical Czochralski puller are shown which will 



steel 




Fig. 1. Geometr^ of a typical LEC Czochralski puller and 
associate global finite element mesh. 
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serve as-a basis for illustrating our method in later 
sections. We will assume at the outset that the 
geometry is axisymmetric. It may therefore be necess- 
ary to distribute in the azimuthal direction some 
parts of the furnace which are not symmetric, such as 
the wiring of the heater. In Fig. 1, the heater H radiates 
power towards the crucible A and the insulating wall I. 
The crucible A rqa^tflins the metal above the melting 
point ; some heat is evacuated through the pedestal E. 
The crystal C is pulled from the melt B with a release 
of latent heat of fusion at the interface, while heat is 
transferred by conduction from B to C. The surfaces 
of the melt and of the crystal radiate heat inside the 
enclosure G towards the wall F which is cooled by 
water circulation. For the case of LEC growth, the 
melt is covered by a boric oxide.layer D. The transfer 
of heat through D is a combination of conduction and 
radiation, since the . oxide layer is semi-transparent 
to infra-red radiation. The semi-transparent property 
"'J will be analysed in Section 3. 

Our main interest is to calculate the shape of the 
liquid-solid interface and the temperature dis- 
tribution within the crystal at various stages of its 
growth. It might therefore be sufficient to limit the 
domain of our calculations to the melt B and the 
crystal C. However, the temperature distribution in 
these "components depends strongly upon .the bound- 
ary conditions on their outer surface, which are a 
priori unknown. Ah accurate calculation of the tem- 
perature field will thus require a global calculation of 
the heat transfer throughout the whole furnace, on 
the basis of its geometry and its material properties, 
and a limited number of control parameters. 

Our method is based on an idealized procedure where 
the crystal is grown with a constant radius selected at 
the outset. This would require an ideal control of 
the growth process ; we generally calculate the power 
input while the growth rate is imposed/ We also 
assume that the whole furnace is in a "quasi-steady 
£ state [4], where one does not take into account the 
0 growth history of the crystal. More precisely, the time 
dependence of the geometry is not taken into account ; 
however, for a fixed situation at a given time, the 
release of latent heat of fusion at the interface is pro- 
portional to the imposed growth rate: The quasi- 
steady hypothesis is valid since the time constants 
corresponding to geometrical modifications are much 
larger than the ones associated with heat transfer 
within the furnace. 

In order to develop our mathematical model, we 
wish to separately consider four types of media within 
the furnace: - ~ 

(i) radiative enclosures -which connect the various 
liquid and solid constituents ; 

(ii) solid domains such as the heater, the crucible, 
the pedestal and other constituents where heat is 
transferred by conduction ; 

(iii) outer thin walls cooled by external convection ; 

(iv) the crystal and the melt together with their 
interface. . - 



. Briefly, the global method consists of separately 
analysing the heat transfer within each of these con- 
stituents and of imposing continuity of temperature 
and:heat .flux on their interfaces. Figure 1 shows a 
typical, finite element .mesh covering . the whole 
furnace. Each separate object (e.g. A, B, C, . . .) con- 
stitutes a macro-element. The analysis will make use 
of several sorts of macro-elements, which are defined 

below. - • * 

A radiative macro-element is an enclosure bounded 
by one or. more macro-elements of another type, 
where heat transfer between walls is essentially radi- 
ative (it is also possible to include convective transfer 
due to the motion of the ambient gas). Temperatures 
and heat fluxes are only evaluated on the walls of the 
enclosure. Since each wall .of an enclosure may see the 
other , walls, , a radiative enclosure cannot be* sub- 
divided into smaller entities. ■ . - \ ;,, : 

Atwchdimensional. macro-element is a liquid or solid 
component where heat transfer occurs by conduction 
(and possibly by convection). It is always allowed to 
partition a conductive (but not. convective) macro- 
element into smaller , entities without ;- modifying the 
transfer properties of the system. .... 

A one-dimensional macro-element is a specific entity 
designed for components which have a low thickness 
as compared to their other dimensions. Typically, the 
walls F of the enclosure G may be well described as a 
thin shell on which a relationship exists between the . 
outgoing flux and the local temperature. ' 
' Parts of the boundary of a macro-element 'may 
either see the outer world, or be interfaced with 
another macro-element. The set of interfaces con- 
necting rnacro-elements constitutes the skeleton of the 
global domain. In the following sections, we will show 
that the global heat transfer calculation consists of 
reducing it to an evaluation of nodal temperatures on 
the skeleton of the furnace. The cost of the calculation 
is essentially related to the number of nodes found on 
the skeleton. 

We will now separately consider each type of 
macro-element while their assembly will be analysed 
in Section 6. 

3. RADIATIVE ENCLOSURES 

In this section, we shall focus on the modelling of 
radiative heat exchange, which plays a major role in 
Czochralski growth because natural emission in the 
infra-red part of the spectrum is important at the 
melting temperature of common semiconductors. We 
will follow the general theory of Siegel and Howell 
[19], although simplifications are required for the 
numerical simulation. Our main assumption is that 
radiation is only - difruse- -This hypothesis , is sat- 
isfactory even 'when specular reflection- occurs,, since 
true reflecting surfaces are generally non-smooth. 
Moreover, we shall consider that emission, absorption 
and reflection of radiating waves only occur at the 
surfaces of the bodies (and not within the bodies them- 
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selves). Further hypotheses will be introduced at a 
later stage. ..-.*- 

3.1. Radiative exchanges in a fixed range of wavelengths 
Let us first establish the integral equation governing 
radiative exchange'in a given enclosure (delimiting a 
volume V) within a nxed r "range A == [A u X 2 ] of wave- 
lengths. We assume that all surfaces are opaque -and 
all material properties are temperature and wave- 
length independent in this' fixed range". In further sub- 
sections, : we shall extend the model " to 'more general 
situations, including the • case of "semi-transparent 
materials. \r. 

The radiation - equation - represents a -balance 
between heat emission, 'absorption and 'reflection. 
Except at the absolute zero, all bodies emit energy in 
the form of electromagnetic 1 waves. In the case "of a 
black body (or a perfect absorber), natural emission 
obeys Planck's law, which has "been determined on the 
basis of "thermodynamic considerations. Let e&iT) 
denote tie spectral emissive" power of the black body, 
i:e." the density of total power emitted by the body 
per unit wavelength at absolute temperature * T and 
wavelength L Planck's law is written as follows': : 



27TC,.. 



(1) 



where C L = 0/595448 kIO 8 W )an 4 m-*..and; C 2 = 
14 388 ^m K are absolute constants! Curves of ^(7) 
as a function of A are given in Fig. 2 for various values 

of r. . v j / : ~ v. ":* ) . 

The total power q^(f) emit^'by the black body 
per unit' area over.. r the" range A of wavelengths is 
obtained by integration of equation (1)' over -this 
range ; one pbtains : . " ' . \~V. 
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Fig. 2. Spectral emissive power of a black body for several 
. . i. * , " temperatures. . ....... .. 
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Rg. S. Fractional black body emissive power in range 0 to 

i , .,jk»P).r j . e M (T)(U-'.: f 

•-• . •■' - . =y K (T)cT\ -..1:1 (2) 

where &'= 5.6696 x'lb -8 W m~ 2 K^ 4 is Boltaiianh's 
constant and y A is defined by the difference 

with _ ■' •• ;' 

15C5d»j i 



i7 5 (e c ^-l); 



(4) 



The curve i^Q is given in Fig. 3. It is clear that y A (T) 
is the ratio between the powers ■. emitted by the, black 
body in the range A and in the -whole. spectrum ]0, cp{, 
respectively, and so lies between O.and 1. Moreover, 
we may observe that y A depends little upon T as long 
as. ■■• ; .. -.T v 

0.2xi0 4 /nnK k XiT<A 2 T< 10 4 urn lL (5) 

In what follows, we shall assume that one of the fol- 
lowing situations occurs.v (i) (5) holds;, (ii) A is the 
whole spectrum (y A = *)> (&) the ^temperature T can 
be estimated everywhere on the enclosure. Thus, in all 
three cases, we assume that y A is a known function of 
the position on the surface of the enclosure. .-' 

Natural emission from a real body is lower than 
that from a black body; As we consider diffuse radi- 
ation and. assume that material properties are tem- 
perature and frequency independent in the range A, 
the total power emitted by the .body per. unit sur- 
face in this range is given by the expression 



• ?ca(x) = £ A (x)? bA (!T(x)) 



(6) 



where £ A (x) is the surface emissivity at the generic 
point x on the enclosure and takes values in the range 
]<U].. . • •* 

Any* element of area is also permanently bom- 
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baxded by waves emitted (or reflected) from other 
points on the same enclosure. The incident heat flux 
$iA(x) can either penetrate the body or be reflected. 
Since we assume the surfaces to be opaque in the range 
A, the non-reflected waves are totally absorbed and 
transformed into internal energy in the body. Kirch- 
hofTs law (which is a . consequence of thermc- 
dynamic equilibrium in isothermal enclosures) states 
that the. absorptivity of any surface, i.e. the ratio 
between the absorbed and incident heat fluxes, must 
exactly equal the surface emissivity. Hence, the 
absorbed and reflected heat fluxes in the range A, 
q^ix) and # rA (x), are given respectively by 

5aA(x) = e A (x)5iA(x:y- (7) 

and.. • . .... . r j.\ . 

g rA (x) = (l-e A (x))g iA : (x), (8) 

As the total outgoing heat flux (or radiosity) ^(x) is 
the sum of the emitted and reflected fluxes at x in the 
range A, we also have 

*«a(x) = q*(x)+q*(x) = £ A (x)y A (x)<r:T 4 (x) 

+ (l-e A (x)) 5iA (x). (9) 

On the other hand, the total incident heat flux q iA (x) 
in the range A is the sum of the contributions of the 
outgoing fluxes. # 0 a( x *) frpm all other points, on the 
enclosure. Let dS and dS* be . infinitesimal areas at 
points x and x*. The fraction, of =the incident flux on 
dS which leaves dS* is .calculated by the product . - 

■ ; . d**M =X(x, x^ oA (x*)<lS < * m (10) 

where K(x t x*) is the surface view factor between x 
and x*. Lambert's cosine law for diffuse radiation 
states that, whenever dS and.^cLS* see each other, 
K(x, x*) is given by the formula 

* (X - X } ~ «[(x*-x) • (x»-x)] 2 - (U) 

where n and n* are the unit 'normals to dS and dS* 
(Fig. 4). On the other hand," when dS and dS* do not 
see each other, Jv(x,x*) vanishes lt 

JST(x,x*) =0.= (12) 

From the knowledge of K(x, x*),' it is easy to integrate 
equation (10) over the surface d V of the, enclosure 




?u(x) 



=1, 



K(x,x*) goA (x)dS*. 



(13) 



Let # A (x) represent the net heat flux at x in the range 
A, assumed to be positive when the corresponding 
energy enters the enclosure.- Clearly, q A (x) is the 
difference between the outgoing and incident heat 
fluxes at x 



^ A (x) = ^ oA (x)-? iA (x). 



(14) 



Calculating ^(x) and ^(x) from equations (9) and 
(14) yields 



9iA(x) = y A (x)<7r 4 (x)- — )-r q A (x) 
«oa(x) = y A (x)ar«(x)-^^^(x). 



(15) 



Thus, introducing equations (15) in equation (13), . 
one finds that the surface temperatures T(x) and 
heat fluxes q A (x) on dV are related by. the integral 
relationship 



g A (x) r 

£a (x) Jfeav 

= ^(x)trr<(x)-£ 



^Cx, x ^^^«a(x*)cLS* 



In order to . obtain a . manageable system of equa- 
tions, equation (16) must first be transformed ..for 
taking into account the axisymmetry of the furnace. 
To this end, let r, 8 and z denote the. cylindrical 
coordinates of x and pipV) stand for the intersection 
of the enclosure surface dV and the mid-halfplane 
$ — 0 (a schematic enclosure is represented in Fig. 5). 
Further, let s be the curvilinear abscissa on p(dV). 
For any x on p(dV) 9 we rewrite equation (16) in the 
form • ••■ • • ■■ i"'.,. i- 



..= y A (x) ff r 4 (x)-| 
J* 



r'^ c (x,x')y A (x0ar 4 (x0oV 
(17) 



where the axisym metric view factor IQfx, x^ is defined 
by the integral 



*c(x : 



X(x,x*)d0* 



(1.8) 



with 



Fig. 4. Surface view factors between infinitesimal areas. 



r x=(r,0,z) 6j p(SF) 

\*=(r',G,z')ep(dV) 

lx* = (r'cos 0*,'r< sin Q*,2T) edV. 



(19) 
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Fig. 5. Schematic enclosure 'section' for an axisymmetric 
furnace (m* = 10). 



Note that, throughout this section, we always rep- 
resent vectors (including position vectors) in Car- 
tesian components; -as: in equations (19); The cal- 
culation of i^x-) from equations (18), (I I) and 
(12) is very complex, since the viewed and hidden 
parts of the enclosure must be taken into account. 
The complete procedure is described in the next 
subsection. 

3 2. Axisymmetric calculation of the shape factor in 
the enclosure 

We wiO assume that the enclosure boundary dVis 
generated by the rotation around the r-axis of 

X J^^ ^ deS Sj * fonDin S the enclosure section 
P(oV) (Fig. 5). Each side is delimited by the edees x, 
and x >+1 , with 5 J 

X J = (r Jt Q 9 Zj). -(20) 
For any pair (x,x-) belonging to the section of the 
enclosure, we need to calculate the axisymmetric view 
factor *;(x, x0, as defined by equations (18) and (19) 
For convenience, we shall assume that, for a fixed pair 
(x,x ), x* is a current point of azimuthal coordinate 
8 related to x' by equations (19). Moreover, n and n' 
will denote the unit normals to the enclosure at points 
x and x', with 



" • > n* = (cos <j>' cos 6*; cbs.^' sin 6* 9 sin 00- (22) 

Whenever x and x* see each 'other, the' view factor 
K(x, x*) _ does not vanish and is given by equation 
(11). Introducing ^ .equanons''(19),7 : (2ly and (22) in 
equation (11), we can express the dependence' of 
AT(x,x*) with respect. tq 0*. We thus rewrite K(x, x*) 
in the form of a function K(8*) ' 

" j^x*) = £(**) = (^^8^^ cos 8-) 
■;• :< • (c+acosfl*)* " • 

• ' : -.J. , , (23) 

where a, Z>, a', and b" only depend on r, z, r'i 

r 7 , # and (which are fixed) but not on 0*. These 
coefficients will not be developed in detail. Integrating 
equation (23) with respect to 8* yields a primitive • 
function 1(8) of K(8*), 



1(8) 



f s 



= A8+B tan- 
sin 8 



■urn—) 



. a+b cos 8 



(24) 



{, 



n : 



' (cos 0, sin <j>) 
(cos <j>\ 0, sin 4> / ). 



(21) 



The unit normal n* at the current point x* is thus 
given (in Cartesian components) by 



where coefficients A 9 B and C are suitable functions 
ofa 9 b t 9 a J t b% a" and b" [12]: * -n : \ 

- -Consider now the integral (18). Let © denote the 
range of values of tf* for which AT(x,x*) does not 
vanish. In order to calculate the axisymmetric view 
factor A;(x,x0 we only need to characterize 0 as a 
set of Z. intervals 0 M/ [ . : . 

' ' 4 * • -v.: . 

U^mi^M/L O<0 m/ < B^^ Ti (25) 

since equations (18), (23) and (24) then provide the 
expected result : 

.: .:; , ^ x ? =^;g ; w^-7(^)). • (2 6) 

It is clear that © depends on the viewed and hidden 
parts of the enclosure. We will therefore analyse this 
last problem in detail. 

Generally, we thus wish to determine whether any 
given line x*x is intercepted of not by the enclosure 
boundary. To this end, we define the circular pro- 
jection p(y) of any point y of the line x*x as being the 
point on the mid-halfplane 8 = 0 having the same 
radial and axial coordinates as y. It is easy to prove 
that any projected, line />(x*x), i.e. the set of circular 
projections of the points forming the line x*x, is a 
segment of a hyperbola (Figs. 6 and 7). We'may 
consider a straight line which is rigidly linked to the 
axis of symmetry and which is initially identical to 
x*x; letting this segment rotate around the axis of 
symmetry will generate a one sheet hyperboloid which 
intersects the mid-halfplane 8 = (f through a hyper- 
bola-passing through x and x' = p(x*). In particular 
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(27) 



is projected on the pair of segments x'x a and x a x, x a 
being the intersection of x' x with the z-axis . 



(28) 



Fig. 6. Projected lines passing through.^ and x for a simply 
connected enclosure section & = 0- The family is indexed 
by the azimuthal coordinate 9* of x*, with x' />(x*). 

the line x^x obviously forms its own projection. More- 
over, the line x'x, where x' is the symmetric of x' 
with respect to the r-axis, as defined (in Cartesian 
components) by - 



Hence, letting 0* (and thus x*) vary in equations (19), 
the set of projected lines p(x*x) will form a bundle of 
hyperbolas passing through x and x' and contained 
in the triangle (x,x\xj. The angle 0* varies between 
0 (if x* ■■= xO and n (if x* = &0 and will be considered 
as the parameter of this set. • 

Besides the fact that they provide a tool for under- 
standing radiative exchange in axisymmetric furnaces, 
the importance of projected lines lies in the following 
observation: for any'line x*x, interception occurs if, 
and only if, the projected line' p(x**Y crossesjhe 
enclosure section p{dV). Equivalent^, a line x«*jls 
not intercepted if, and only if, its projection p(x*x) 
does not cross any side Sj ofptfV). This property will 
allow us to calculate the axisyrnmetric' view factor 
Kcfr x% For any side Sj ofp(dV) (Fig. 8), we define 
the interval [6 mJ , 6 MJ ] as the range of values of the 
azimuthal coordinate 0* of x* such that the projected 
line p(x*x) is intercepted by Sj 



0*e[0 ro ,,^^*(x,x*) = O. 



(29) 



Note that these intervals may be empty and may also 
overlap. The range © of values for which K(x,x*) 
does not vanish is thus given by. • 





Fig.' 7. Same legend as in Fig. 6 for a multi-connected 
enclosure section (£, = 2). 



Fig. 8. Interception of a family of projected lines by a side 
(Sa) of the enclosure section. 
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®^OM\jj i [S aJ J w] . (30) 

Since knowledge of © provides *V(x, xQ from equa- 
tions (26) and (25), it is clear that our final problem 
lies in determining the intervals [B mJ , 6 MJ ). These can 
easily be calculated by applying the theorem which 
follows. 

Theorem. Let Sj be any side of the enclosure section 
S «tremities x, : and x y+ , given by equation 
(20). Let also x be a fixed point on p(dK) and x* be 
on the enclosure dV, a point of unknown azimuthal 
coordinate 6* and fixed circular projection x' (as given 
by equations (19)). Two cases are considered (Fig. 8) : 

(i) . if the projected hue p$*x) passes through 1 x- 
then 6* can be obtained by solving the equation 

(ii) if X^x) is . tangent to -.^i, then 6* can be 
Obtained by solving the.system 

o* = e+e:. • 



P - £2k± 1 " ± ( r ' z J+ i-r J+l z,) 
cos0' = z 'frV+ 1 -O) + (r y z /+ 1 - r<4 , , Zf ) 



(32) 



When they occur, indeterminacies must be removed 
by application of L'Hospital's rule. In both cases, the 
searched projected line exists if, and only if, equations 
(Ji; or (32) have a real solution. Moreover in the 
second case, the coordinates of the contact point K- 
when it exists, are given by the formulae 



■ 1 



r c = 



rr'sin(fl-)-fl') 
rsinfl+r'sinfl' 
zr' sin d'+z'r sing 

r sin 0+r' sin g' " 



(33) 



The proof of this theorem is solely based on geo- 
metrical considerations and will be omitted. 

r4 Le tf US , n £T retum 10 the calculation of the intervals 
\Pmj, V MJ l The previous discussion shows that S , and 
By must be chosen as the lowest and highest value of 
f , respectively, between the following: 

K J! **r = v 9r , (and) °* = *> * S J "to?** one or 
both ofjhe degenerate projected lines ?x" and 
x'x,vx,x; • 

(ii) the values of > obtained by applying the pre- 

both of these are located inside the triangle (x, x. xV 

(iii) the value of 8* obtained by applying the'pre- 
vious theorem (second case), if the contact point L is 
located between x, and.x, + 1 , and inside the triangle 



3.3. Discretization of the radiative integral equation 

Any radiating enclosure of the furnace is discretized 
by means of a set of n v one-dimensional finite' 
elements. The approximation of temperatures and 
heat fluxes is everywhere continuous on the enclosure 
section p(dV), except atits edges (where analytical 
heat fluxes are normally discontinuous). Using quad- 
ratic shape functions, we write in equations (17) (for 
wavelengths in the range A) 



(34) 



'?a(x) = 2>^,(x) . 

We thus approximate 

= (l^(x)J by ZTfMx). 

Numerical experiments have shown 1 the validity of 
such an approximation which considerably simplifies 
our calculations. 

We discretize the integral equation (17) on the basis 
of a Galerkin formulation, this method has been 
preferred to a procedure'of nodal collocation because 
m the latter case, important" numerical problems arise' 
from the calculation of the axisyrrunetric view factor 
£ c (x,x / ) at the edges of >(5V). thus, discrete equa- 
tions are obtained by multiplying both sides of equa- 
tion (.17). by (pfrfa) ds) and integrating the result over 

.*■"- . . yr . .... . 



c*p) ... e A (xO 



x^(xOaVd?Jcrr;. 



(35) 



After suitable matrix operations, the system (35) 
may then be written in the form 



(36) 



where q A and T 4 respectively denote the column vec- 
tors of nodal heat fluxes (in the range A) and fourth 
power of temperatures, and Fj, is now a known matrix. 
It is clear that either A extends over the whole spec- 
trum ]0, oo [, and g A (x) is the total heat flux q(x) enter- 
ing the enclosure, or several ranges of wavelengths 
Ai,A 2 , ... must be considered. .In the former case, 
dropping the index A in equation (36), we write 
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. q = £T, (37) 
In the latter case, we have 

= '(^ 1 +£ 3 +--0T 4 (38) 

and thus the total heat flux depends on the fourth 
power of temperature as before. Equation (37) is the 
general discretized relationship governing radiative 
exchange in an enclosure. The case of semi-trans- 
parent materials also reduces to this form and will be 
investigated in the next subsection. 

Numerical integration of equation (35) with respect 
to ds is performed with a three node Gaussian quad- 
rature on every element, while 20 nodes are used for 
the" ; integration with respect to oV. This choice is 
justified by the steep variations of .K^xO observed 
in the vicinity of the edges of p(dV}. However, even 
^ with such an accurate method, poor results have been 
' obtained in some circumstances. It is possible to check 
the accuracy, of the integration with the help of the 
identity . . 

In the worst cases, the left-hand side was equal to 
0.3 •■• - 0.2. We have found that the origin of this prob- 
lem lies, in the fact that the axisymmetric view factor 
is not a regular function, whereas' this is required -for 
applying Gaussian iquadratufe: . We 'have therefore 
improved the algorithm in the following way : as equa- 
tion (35) is integrated with respect to both ds and dY, 
the integration with respect to ds' must be performed 
for any given integration node x (recall that three such 
nodes are located.: on any element). For any fixed x, 
we divide the enclosure section p(d V) at any point x' (1 ^ 
x^,-. t . where ^(x, x 7 ) is a discontinuous function.of 
S* , or has a discontinuous, derivative with an infinite 
slope on one or both sides of the discontinuity. For 
Qy any element where : such a division point is located, a 
20 node Gaussian mtegration (with respect to dV) 
is then performed separately on both parts of the 
element. . v.- • . 

It is possible to prove that K c (x,x') is not regular 
fa 1 : previous sense) .w hene ver one of the 
dejgenerate lines xx (0* = 0) or xX u x tt x (0* = it) 
passes through one of the edges of the enclosure 
section. These cases can be easily detected (Figl 9). 
On this basis, the improved algorithni can provide 
accurate results as long; as the identity (39) is used 
to check the validity of the mesh of the enclosure: 
generally, when two elements . are adjacent , to a -re- 
entrant edge ofp(dV), their respective sizes must be 
of the same order of magnitude. ■ : . r : 

We conclude this subsection by noting that Galer- 
kin's method is globally conservative, which means 
that the total heat £ux balance would vanish on the 
enclosure if integrations were performed exactly. As 
a matter of fact, the identity '■ 




Fig. 9. Degenerate projected. lines passing through one of 
... the edges of the enclosure section p{dV); 



• E ^(x) = 1 - (40) 

can be introduced in equation (35) through the sum 
of the latter equations (which are indexed by /)• 
Reversing : the ' order of integrations, taking into 
accounjt the symmetry of the kernel K c (x,x r ), and 
using the identity (39), one ; can then obtain 

.i'H:'^^^^" 0 - (41) 

or equivalently, ; using equations .(34) 
*.' \ :^ A (x*)dS'*- 0. - (42) 

3 .4. Semi-transparent materials, . ..... - : , , 

. The rboric.. oxide layer used, m the. Czochzalski 
growth of gallium arsenide for avoiding arsenic evap- 
oration is not opaque and thus one of the hypotheses 
specified in the first subsection does not,.hold. In a 
semi-transparent medium, absorption does not take 
place solely at the surface of the body; Energy is partly 
transmitted through the body, and also. absorbed and 
emitted within the body itself. 

Let us consider a monochromatic ray of wavelength 
A, with an initial intensity /J , which crosses (without 
reflection) a material of thickness h: It is attenuated 
by absorption and dispersion. Neglecting the latter, 
Bouguer's law gives a relationship between the out- 
going intensity I h x and the absorption coefiicient a x • . 

. 0 . . .. A = It exp . (43) 
The penetration depth L x is defined as. . . 
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(3) 



(3) 



(1) 



, (2) 



'sem i- transparent J 
'< layer 



(3) 



(1) 



(1) 




+ 




Fig. 10. Application of the band-energy method for calculating heat exchange in a semi-transparent 

• material. 



Comparing this depth .with the thickness h of the 
material, one obtains the spectral optical thickness z x 



(45) 



fe(x) = ^ Al (x) on(l) .. 

\ ?(*) « Sa 2 (x) :r ; on (2) ^ '(46) 
;; ,U(x) ^ ^XxK^tx) "on (3) v " ! " : 
It is therefore possible to obtain a general law of the 
form of equation (37) over the whole boundary by 
superposing (as in equation (38)) the contributions 
of each enclosure. The resulting enclosure may be 
considered as a single, radiative macro-element. 

4. TREATMENT OF TWO- AND ONE- 
DIMENSIONAL MACRO-ELEMENTS 

. In this section, we wish to study macro-elements 
where heat is transferred by conduction, such as the 



heater and the crucible, with the exception of the melt 
and the crystal which will be considered separately in 
Section 5. .In such macro-elements, the. temperature 
field is governed by the equation - 



When z x is large compared to one," the material is 
transparent for the wavelength X; when z x is sm all 
compared to one, the material is opaque for the wave- 
length A; Our model assumes that the ranges of wave- 
lengths for which z x is of the same order as one are 
negligible; using the band-energy method, we thus 
consider that the boric oxide layer is transparent for 
some ranges of wavelengths (z x > 1) and opaque for 
others (z x < 1). '[ ' 

On this basis, it is easy to calculate radiative' heat 
exchange in a semi-transparent medium. Let us con- 
sider a simplified example, represented in Fig. 10,' 
where we have a semi-transparent layer. We define 
two different enclosures, where radiative exchanges 
take place with different wavelengths'. Enclosure I cor- 
responds to waves On the range Ai) for which the 
layer is transparent, while enclosure II corresponds to 
waves Cm the range A 2 ) for which it is opaque. Within 
each enclosure, radiative exchange is governed by the 
integral equation (17), which is approximated by 
equation (36). We separate the boundaries of the 
domain into three groups, denoted by "(1), (2), (3) 
(see Fig. 10) : boundaries (3) are 5 common to both 
enclosures, while boundaries (1) and- (2) belong to 
enclosures I and II only. The total radiating heat flux 
q(x) is thus given by 



" V • \k(x) V7Xx)]^ r(x) = 0 



(47) 



where 7c is " the i thermal conductivity which may depend 
upon the position x and r is 'a heat source' per unit 
volume which will in general depend upon X. For 
simplicity, we assume that k is temperature inde- 
pendent in order to preserve the linearity of the 
system. We will first discuss the case of two-dimen- 
sional macro-elements: -'J ... .- r .. ; 

The macro-element is covered by- a* mesh of finite 
elements, which in, the present paper are six-node -tri- 
angles or nine^node quadrilaterals. The ;apprpximate 
temperature field is represented by the sum 



\ (48) 



where the 7ys are nodal temperatures forming a vector 
T and the 0/s are the corresponding global shape 
functions. Let O. be the axisymmetric domain cor- 
responding to the macro-element; applying Galerkin's 
method for discretizing equation (47) , one obtains 



*f[y : fcY'D+r]dp : £iO 



.(45) 



or, with an integration by parts 



£ ktyfr -'VT) dQ £ £ ^r dn-^ frq oT (50) 



where X is the boundary, of Q and q t is the outgoing 
heat fiux, as given by , 



(51) 



n being the outer normal on £. The boundary 2 is the 
union of two parts, £ + and S° (Fig. 11) ; 2 + belongs 
to the skeleton of the furnace defined in Section 2, 
while X 0 is the part of the outer boundary where the 
temperature is fixed or is related to the heat flux den- 
sity q by a law such as 



q = h(T-T 0 ) 



(52) 
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Fig 11 Boundary partition for the calculation of heat trans- 
fer in a conductive macro-element 

where A is a convection coefficient (W m" 2 K*jj) and." 
Tp is the coolant temperature. , : 

Let Tt denote the vector of nodal temperatures on 
£ + and let T* denote the vector of. all other nodal 
temperatures associated with the macro-element ; we 
write 

• J $ -^yT^jr+.T*].' r-'.,:- ... . " ' (53) 

It is then easy to. rewrite equation (50) in me following 
form;' ; ■ " "' 

iSr T T*+M + T + =E + (»0-Q + ' (54) 
where the folio wing, symbols have been introduced; . . 

~ " .{£ ***** * ^1 " ' • ' " "• 

Q+ -{J r+ ^;^} : - (55) 

The superscripts * and + refer to the partition of nodal 
temperatures, as denned in equation (53). The symbol 
Win equations (54) represents the total power input 
in the macro-element which, at -this : stage, may.be 
unknown. ; . ... ; - 

In order to link at.a later stage the two-dimensional . 
macro-element with the adjacent radiative enclosures 



or with any other macro-element, we need a rep- 
resentation of the heat flux density q along 2 + . The 
latter is assumed to be " • ■ ■ 



(56) 



in terms of nodal flux densities and shape functions.. 
The number of nodal flux densities along 2 + ; may be 
larger than the number of nodal temperatures. Indeed, 
if E + has corners, we will associate two nodal heat 
flux densities at these corners, one for each adjacent 
section of the skeleton (Fig. 12)';' we will show in 
Section 6 how 'to handle such multiple values in the 
global calculation on the skeleton. Let denote the 
vector of nodal heat flux densities on S + j its dimeh- . 
sion may thus be larger man that 6f T+. (^mbining 
equation" (56) and the sixth /equation (55), one obtains 
a relationship of the form 1 



(57) 



which is "introduced m equation (54) . 

At this stage, it is possible to perform a jr/aric con- 
densation of the macro-element or, in otter words, to 
eliminate T* from the system' (54). Formally, if one" 
defines the matrix ; . i 



and the vector.. 



(58) 
(50)- 



,(6Q) 



one obtains the system , i • 

Such a static condensation' is/easy to perform when' 
the frontal elimination method is 'used. We have thus - 
shown that it is possible to reduce a heat conducting 
macro-element to a linear relationship between the 
nodal temperatures and the nodal heat fluxes on that 
part of its boundary which belongs to the skeleton. 
The system will be closed when the various macro- 
elements are connected by means of the radiative 

enclosures. ,' 

We have seen in Section 2 that, for some parts 1 of 
the furnace, it is interesting to use one-dimensional 




Fjg. 12. Nodal flux densities and temperatures at the corners 
■i . of the boundary. 
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radiative 
enclosure 



q - h (T-T 0 ) 



Fig. 13. Example of a one-dimensional macro-element on 
the upper dome of a Czochralski furnace. 



macro-elements. An example is given in Fig. 13 : the 
upper radiative enclosure is adjacent to a thin shell 
(representing the dome) which is part of the skeleton. 
The loss of heat alo^g the outer side of such a macro- 
element is again given by a relationship of the type 
of equation (52) ; since me niacro-element does : not" 
accumulate heat, we may thus write "" 1 ' ■ 



ff + +/t(7?-r 0 )«o 



(61) 



where q+ denotes the heat flux density towards the 
skeleton. Using Galerkin's method : with a : rep- 
resentation similar to equations (48) and (56), it is 
easy to obtain an equation of the form* . 



- (62) 



forthe one-dimensional macro-element. 

Let us close this section by mentioning a particular 
use of one-dimensional macro-elements which is 
found useful .for dividing a radiative, enclosure in two 
entities with negligible interaction between them; an 
example is given in Fig. 14. Two adjacent one-dimen- 
sional macro-elements are then used, with the con- 
straint that nodal temperatures at adjacent nodes are 
identical while the sum , of the corresponding nodal 
heat flux densities vanishes. 



double ID 
macro-element 



\ 



0 



5. THE CRYSTAL-MELT MACRO- ELEMENT 

In this section, we wish to study the domain formed 
by the crystal and the melt It will be found important 
to associate these two components within a single 
macro-element, because they are separated by the 
liquid-solid interface which is a priori unknown. The 
interface originates at the point of intersection of the 
solid, liquid and gas (or oxide in LEC growth) 
domains, or tri-junction, where the temperature must 
be equal to the melting temperature of the crystal 

Locating the interface is more difficult with the 
Czochralski growth (see, e.g. ref. [34]) than with the 
Bridgman process (see, e.g. ref. [35]). In the first case, 
an appropriate combination of pulling rate and power 
input must be calculated in order to grow a crystal 
with a constant predefined radius ; we will show that 
the optimal technique consists of imposing the crystal 
diameter and the growth rate," and of getting the heater 
input power as a result In the second case, the inter- 
face is simply defined by the melting isotherm, and 
• the tri-juhction is thus an additional degree of free- 
dom; this means that the pulling rate and the input 
power may "be' imposed, while the interface location 
results from the calculation^ 

5.1. Basic equations ~ ;r 

A cross-section of a schematic crystal-melt macro- 
element is shown in Fig. 15. The axisymmetric domain 
CI is the union of a solid region H s arid a liquid region 
Q L , separated by an interface T. Our mathematical 
model relies on a set of hypotheses which we repeat 
for the sake of clarity : 

(i) the quasi-steady state hypothesis allows us to 
neglect the partial time .derivatives, in . the governing, 
equations; , 

(ii) the length and the radius of the crystal are 
imposed but the location of the interface is unknown ; 

(in) the shape of the free surface of the melt is 
imposed at the outset of the global calculation; in 
particular, the shape of the meniscus linking 1 the crys- 
tal to the surface of the melt is calculated separately 
(see Section 5.4) ; ***: r . ..." • ' ... 

(iv) additionally, we will neglect the change of den- 
sity from the liquid to the solid domain. 

When the motion of the melt is taken into account, 
heat transfer in the liquid domain is governed by the 
following equation : 



Fig. 14. Use of a twin one-dimensional macro-element for 
subdividing a radiative enclosure in two parts. 



pcv-Vr-V-foVT) = 0 



(63) 



where v is the velocity field in the melt, p the-density, 
c the heat capacity and k L the thermal conductivity. 
In the solid domain, we have • ■ 



v--(* s vr) = 0 



(64) 



where k s is the thermal conductivity in the crystal : In 
view of the low growth rate of the crystal, vertical heat 
convection can be neglected by comparison with the 
diffusive term. ... 
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Fig. 15. Cross-section of the crystal-melt macro-element 



In order to evaluate the release of latent heat of 
fusion along the interface T, let us denote by v n (s) 
the normal growth velocity of the crystal along tie 
interface as a function of an arc length s, with respect 
to a frame attached to the seed; thus, i> n is positive 
during growth. The amount of solidification per unit 
time and area along T is the product pv U9 and the 
corresponding release of heat is given by - • ■ 

W r (s) = P AHrV n (s) (65) 

where AH { is the specific latent heat of fusion of the 
crystal. Under quasi-steady , assumption, the normal 
velocity v n is related to . the pulling rate of the 
crystal by the relationship , , -.■ 

= -»p»jns v e/ ^ (66) 

where e z is a unit vector .pointing upwards and n s is 
the outer unit normal to fig along T. 

Denoting by g L and g s the local outgoing heat flux 
densities along T, we require that .the energy balance 
be preserved along the interface and thus 

^rW+9LW + ft(j) = 0. • ; - ; (67) 

If one neglects supercooling, the location of the 
interface itself is found through its identification with 
the melting point isotherm, i.e. : 



71G0 = T m , along T. 



(68) 



5.2. Spatial discretization 

.•.•Since the location of T is a priori unknown, the 
discretized form of the basic equations must be estab- 
lished on a deformable finite element mesh which fol- 
lows the -evolving form of the -melting isotherm 
throughout the iterations. This method was first intro- 
duced for solidification, problems by Ettourney and 
Brown [36]. Figure 16 shows a typical crystal-melt 
macroelement, where the shape of the meniscus has 
been calculated at the outset and where the location 
of the tri-junction is imposed; As an initial guess, one 
assumes that the interface is fiat. The domains Cl s and 
Cl L are then covered by an array of quadrilateral nine- 
node Lagrangian elements. These nodes are supported 
by a set of vertical lines. Below the free surface of the 



melt, the altitude of the nodes is fixed during the 
iterations ; on the contrary, the altitude of the nodes 
below the interface and within the crystal will vary 
with the position of the interface, in order to maintain 
geometrically well-behaved finite elements through- 
out the iterations. The various techniques for man- 
aging moving meshes automatically have been 
reviewed in ref. [37]. • ■ 7 • ' • 

Let z t denote the altitude of a node i on the 'interface. 
The location of the interface is then fully identified by 
a vector z containing the z/s of all the nodes of T. Let 
k denote a node lying in. the melt below node i; its 
altitude z^ will then be given by 

where c fa - is a proportionality factor which is imposed 
at the outset. The same rule is used for locating the 
nodes in the crystal. ..- . . - . 

The concept of deformable mesh has important 
implications on the finite element discretization of the 
temperature field, which is now represented, in H L as 
well as Qs, by a sum of the form - . > : 



(70) 



where we note that the shape functions depend upon 
the location of T through the position of the nodes. 

Applying Galerkin's method as we did in Section 4 
to the balance equations (63) and (64), one easily 
obtains 



I 



J^rfdS- J r *,FK r dT (71) 



since*. from equation (67), we have 

J r ^L+ffs)dT=-J r ^^ r dr. (72) 



On the outer boundary 2 of £2, the heat flux density 
has a representation similar to equation (56) while, 
along r, the nodal temperatures must satisfy equation 
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' " Fig. 16. Typical moving mesh for the costal-melt macro-element. 



(68). In what follows, we will not take convective 
terms into account ; their impact will be discussed in 
a later paper (see also refs. [15, 17]). 

We will now use the same decomposition as in 
Section 4 and write for the crystal-melt, macro- 
element- 

• T«rrvr+] craj 

where T* is the vector of internal nodes of the macro- 
element and T + is the vector of nodal temperatures 
on its boundary. It is then easy to rewrite equation 
(71) as follows: 5 

M*(z)T* +Mz)T + = E*(z, M . 

^ t (z)T*+M + ( z ) T+ +^ f ^ + = E+ ( z ^pui) < 74) 
where M*, M*, N, B + and q + are denned as in equa- 
tions (55) and (57), and where 



E + (z,u pul ) ~ 



W r dT\ 



■} 



(75) 



We observe that the matrices M*, and N depend 
upon the location z of the interface while E* and E + 
also depend upon the pulling rate. For a given value 
of z and i?^,, it is again easy to perform a static con- 
densation and obtain 

^ + (z)T+-r £ + q + - C + (z, t> pui ) (76) 
where A+(z) is denned as in equation (58) while we 
have 

C* (z, i^) = E + (z, v pvl )-N T M*~ 'E*(z, 0pal ). (77) 
The vector z is unknown at the outset; the cor- 
responding system of equations is 

T r =T t - (78) 



where T r is the vector of nodal temperatures along T, 
and T t is a vector with all components equal to the 
nodal temperature at the. tri-junction. , 

5.3 . Local iterative scheme 

We will show in Section 6 how the calculation of 
the temperature field and the z vector within the crys- 
tal-melt macro-element can be decoupled from the 
calculations on the rest of the furnace. Let us assume 
that either the temperature or the heat flux is specified 
at each node of the boundary of the macro-element 
The non-linear system of equations (74) and (78) is 
then solved by means of Newton-Raphson's method. 
Let T*, T + , q + and z denote the set of unknowns after 
a given iteration ; the corrections 5T*, £T + and Sz are 
then obtained through the solution of 
M*FT* + dz • [V r M*T* + V Z NT+ -V r E*] 

a E*— M*t*-N? + ^ 
tfW + Sz ' [V X ^ T T + - V*E + ] ... .-. • 

$T r = 0. .(79) 
■ * . - ■ ." 
The iterative procedure is interrupted once the mag- 
nitude of the right-hand sides in equations (79) lies 
below a pre-assigned criterion. 

In developing equations (79), we have assumed that 
the pulling rate v p * was known. In Section 6, we will 
discuss an algorithm where r pul is unknown ; it is then 
necessary to add to both left-hand sides of equations 
(79) a term of the form : 

(i) -^puj^E*/^/ . 7 . 

(ii) -5t7 P ui 5E + /^pui5 " 

respectively, while the equation corresponding to the 
new unknown c pul states that the tri-junction is at the 
melting temperature 
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Fig. 17. Fixed wetting angle <f> between the vertical wall of 
the crystal and the surface of the melt 



r t = r m 



(80) 



5.4. Calculation of the meniscus shape 

Before closing this section, it should be recalled 
that, for an accurate calculation of the melting point 
isotherm and of the heat transfer near the tri-junction, 
it is necessary to know the shape of the meniscus 
which links the edge of the crystal and the upper 
surface of the melt 

• An important parameter for calculating the men- 
iscus shape is the angle 0 shown in Fig. 17 between 
the vertical surface of the crystal and the surface of 
the melt. ~ We assume that this wetting angle <j> is a 
fixed material property. The latter has been postulated 
by Bardsley et al. [38] (through arguments of ther- 
modynamic equilibrium) and has been verified exper- 
imentally by Surek and dhalmers [39] for silicon over 
a wide range of growth rates. For the case of a crystal 
obtained by the LEC method, one should also con- 
sider a second meniscus at the boundary between the 
crystal and the upper surface of the boric oxide layer. 
For this second meniscus, we also assume that the 
angle of contact can be imp'osect. 

The locations of the melt and encapsularit menisci 
are governed by the Young-Laplace equation' expres- 
sing a force balance between surface tensiori'and grav- 
ity [40]. The equation is solved at the outset by a 
decoupled finite element scheme * the shapes of the 
menisci are left unaffected during the global iterative 
procedure. 

6. GLOBAL SOLUTION 

We have shown in earlier sections that, for every 
kind of macro-element, it is possible to obtain a system 
of algebraic equations connecting nodal temperatures 
and heat fluxes along the skeleton. To summarize, let 
us assume that the puller may be associated with n T 
radiative enclosures. For each of these, we obtained 
in Section 3 a system of the form 

q<*>«F<*>T<*>* l<Jb<* (81) 

where the full matrix F (k) connects the nodal heat 
fluxes to the fourth power of the nodal temperatures. 
Let us also consider n c conductive macro-elements, 
including the heater as well as the crystal-melt macro- 



element. For each of these, we found in Sections 4 and 
5 a relationship of the form 

^(z)T^-h^q^ = C^(z, W), U;^ 

.(82) 

where again T 01 and q 00 are the nodal temperatures 
and heat fluxes on the skeleton part of the boundary 
of the domain. The dependence of A9} and upon 
z and t? pul only occurs for the crystal-melt macro- 
element, while W only appears in - the • macro- 
element^) where heat is being dissipated.- We note 
that Ay* has a block structure while 2? w 'has a band 
structure. . ■: 

Equations (78), which correspond to the unknowns 
characterizing the position of the interface, must also 
be considered. Moreover, since we have assumed that 
the radius of the crystal is fixed at the outset, the 
pulling rate t> pu j and the input power FK.are dependent 
variables ; however, one or the other is assigned as. part 
of the data for a global simulation. The "equation corre- 
sponding to the remaining unknown is equation (80). 

The system (81), (82), (78) and (80) is highly nonlinear, 
in view of the fourth power of T in equation (82) and 
the dependence of the various matrix operators upon 
z. The need for an iterative procedure is obvious. 
Decoupling is also necessary in view of the size of. the 
system. On the basis of , in equation (82), it is 
possible to recover the temperature field in any macro- 
element by means of adequation of the.form 

^*</)(z)T*^+7^(z)T^ - £* w (z,r puJ , W). (83) 

6. 1 . Skeleton equations 

Before describing the global iterative algorithm, let 
us asume that z, v pxjl and W are known, and let us 
show how it is then possible to calculate the tem- 
perature field on the skeleton, without using equations 
(78) arid (80)! The term 'B&qU> in equation (82) stands 
for generalized nodal heat fluxes Q 00 along the'bound- 
ary of the macro-element (see equation (57)) . Along 
the common border of two ' conductive ' macro- 
elements, such nodal heat fluxes are of opposite sign. 
They can thus be easily eliminated from the global 
system simply by adding the appropriate equations of 
the form of equation (82) along the common borders. 
It is thus possible to obtain a system of the form 

: . . A(z)T+Bq* - C(z, n^ 9 W) . (84) 

which relates the nodal temperatures T s and the radi- 
ative heat fluxes q r along the skeleton. We note that 
the matrix B is rectangular since the nodal heat fluxes 
have been eliminated alon£ the common border of 
adjacent conductive macro-elements. In the same way, 
we write equation (8.1) as follows : 

q r = FT* 4 (85) 

where Fis also rectangular. 

At geometrical corners of radiative enclosures, we 
had been led to identify a different nodal heat flux on 
both sides of the corner, while the associated tern- 
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perature is obviously unique. For handling the system 
(84) and (85), we have found it easier to also define 
two temperatures at geometrical corners in order to 
have the same number of nodal heat fluxes and nodal 
temperatures. The equality of both temperatures at 
the nodes is then imposed as an additional constraint. 
Introducing equation (85) in. equation (84), we obtain 

A(z)T+BFT A = C(z, v paU W). (86) 

Assuming again that z, and W are known, this is 
a non-linear system in terms of the skeleton tem- 
peratures. Jt is solved by means of the modified New- - 
ton's method (also called Newton-Richardson [41]). 
A LU factorization of the Jacobian matrix is executed 
by means of an .algorithm which is well suited for 
hollow matrices with irregular structure [42]. Having 
solved equation (86), we obtain the set of nodal tem- 
peratures along the skeleton. From these, it is then 
possible to calculate the nodal temperatures inside 
macro-elements by means of equation (83), while radi- 
ative heat fluxes can be* retrieved from equation (85). 

6.2. Global algorithm \W imposed 

Let us first assume that, for a given crystal radius, 
the input power Wis imposed while the corresponding 
pulling rate is unknown. For solving the system^ we 
use a decoupled iterative procedure which is described 
as follows. Let z„ and u^. denote the corresponding 
values of z and after iteration n ; typically, we will 
assume that Vanishes while Zq corresponds to a 
horizontal interface^ 

On this basis, we solve equation (86) and obtain a 
set of nodal temperatures T* +1 along the skeleton 
from which we calculate, with the help of equation 
(85), the nodal lieat fluxes q T „ + , . 

Next, we perform a decoupled calculation on the 
crystal-melt macro-element. for calculating z„ +1 and 
i$S With that objective in mind, we solve the system 
of equations (74) and (78) with t? pul as an additional 
unknown ancl equation (80) as an additional con- 
straint; as boundary conditions, we use TJ +1 along 
•the melt-crucible interface, and , along the surface 
of the melt and the wall of the crystal. The new values 
of z„ + , and i£J 1 are then used for. recalculating the 
matrices in equation (86). the iterations are pursued 
until convergence criteria for T, z and are satisfied. 

The above decoupled iterative procedure converges 
in general, but can lead to unexpected and often un- 
realistic results. Indeed, it is found in practice as well 
as in theory that a slight modification of the power 
input can lead to . important variations of the pulling 
rate for keeping the crystal diameter constant Since 
the power input cannot be precisely measured on 
industrial pullers, one can easily obtain unrealistic 
huge values of the pulling rate. 

6.3. Global algorithm : v pu , imposed 

In view of the above considerations, we found it 
much more efiicient to impose the pulling rate, which 



is a control variable in industrial growth, and to cal- 
culate the necessary power input. The decoupled iter- 
ative procedure may be described as follows. Let z n 
and W n denote the corresponding values of z and W 
after iteration 72; we will assume for W 0 a first guess 
based on a rough estimate. We solve equation (86) 
for T and W 9 with equation (80) as an additional 
constraint, i.e. we select the power input such that the 
temperature at the tri-junction is precisely the melting 
temperature. We obtain T„ + 1 and <£ + , and then per- 
form a decoupled calculation on the crystal-melt 
macro-element for calculating z„ + , while r puI is fixed. 
The boundary conditions for the decoupled problem 
are the same as in Section 6.2, except that constraint 
(80) does not apply at this stage. The iterations are 
pursued until convergence criteria for T s , z and Ware 
satisfied. 

^ The stability properties of the algorithm where t? pu] 
is imposed are remarkable, and definitely better than 
those found when W is imposed. In practical appli- 
cations,, a -relative accuracy of l%o is obtained with 
six global iterations; The . solution of equation (86) 
requires some four iterations for the .first resolution 
and only one for the next iterations, while, equations 
(74) and (78) require some four- iterations.: 

7. EVALUATION OF THERMAL STRESSES 

Performing the fully coupled calculations of Section 
6, we are able to obtain an* accurate' temperature field 
and the shape of the tfquid-^solid interface during 
growth. We may thus calculate the stress field without 
any arbitrary simplifying assumptions on the shape 
of the crystal and the temperature, field. A' realistic 
calculation of the stress field would in turn require a 
continuum "model taking Mo account the material 
properties of the crystal near the melting temperature. 
However, in view of the lack of material data at such 
high temperatures for the metals which we have con- 
sidered, we have been led to formulate the following 
hypotheses:. " V '"" 

-(i) the crystal is an isotropic linearly elastic solid 
and the stress field is axisymmetric ; deformations are 
reversible and plasticity and creep are not takeninto 
account; 

(ii) the thermal field is decoupled from the stress 
field; ' ' ; >" 

(iii) the material, properties of the crystal (Young's 
modulus, Poisson's ratio, thermal expansion coefficient) 
are temperature independent ; 

(iy) the crystal is free of surface forces ; ' 
(V) the crystal is stress free under a uniform tem- 
perature field.. • r 

The finite element calculation of thermal stresses 
under such conditions is straightforward (see, e.g. 
ref. [31]): A quantity of interest for evaluating the 
generation of dislocations is the Mises invariant which 
is defined as follows : 

s M -= Ks 1 ^s 2 ) 2 +(s 2 ^s 3 y-h(s 2 -s l y]^ 2 (87) 
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where S u S 2 and S 2 are : the principal stresses. It 
is generally accepted that plastic deformations occur 
where S M reaches a critical value. 

: 8. RESULTS 

In the present section, we wish to demonstrate the 
efficiency of our numerical method through an analy- 
sis of two typical calculations. The first set of results 
corresponds to a Czochralski furnace for the growth 
of germanium crystals; we will show that, for cal- 
culating the growth of a crystal with fixed radius, it is 
much more efficient to impose the pulling rate than 
the power input. Our second example refers to the 
LEC growth of a gallium arsenide crystal ; we will find 
that the semi-transparent properties of the oxide layer 
have a major impact upon various aspects of the 
solution. • 

8. 1 . Germanium growth 

A typical furnace for growing - germanium crystals 
is shown in Fig. 18 ; for the sake of clarity, a different 
scale is used in the radial and axial directions. In the 
present, example, the crystal ^diameter is . 10 cm, its 
length is 20 cm, and the total weight of metal is 30 kg. 
The meniscus near the tri-junction has little impact in 
the present calculation . .and has been neglected. The 
germanium melting temperature is 121 1 K ; the other 
physical properties of the materials are giyen in Table 

It . should be noted that some thermal properties, 
and the emissivity in particular, are highly dependent 
upon the grade and. the surface, state of the material. 
This explains why some values for the same material 
appear to be different for both germanium and gal- 
lium arsenide furnaces. 

The finite element mesh, used for the calculations is 
shown on the .left part of Fig. 19. The global mesh 
contains 965 biquadratic elements, and ,43.17 nodes. 
im The skeleton supports 792 nodes. A typical growth 
[jfr) rate for .the process is .3 cm h" 1 .; the, corresponding 
isotherms are shown on the right part of Fig. 19. 
The crystal is slightly concave at the interface, with a 
deflection of 0.6 mm The necessary power in the 
heater is 24.5 kW. To obtain the solution, we have 
performed six global iterations with adjustment of the 
power input; for each global iteration, we needed 
three inner iterations for locating the interface; the 
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Fig. 18. Draft of the germanium furnace used for the simu- 
lations (radial axis is dilated). 



convergence criterion was 10" 6 K for the temperature 
field and 10~ 6 m for the position of the interface. 

An interesting problem is to determine the 
maximum value of the growth rate for a given process 
stage. In view of the high non-linearity of the problem, 
the pulling fate can only be increased in this case by 
a series of mcrements (the convergence of the iterative 
procedure could not be obtained otherwise). Let us : 
examine in Fijg. : . 20 the isomernis which have been 
obtained for values 5 of t? pu i equal respectively to 4.2, 
4.4 and 4.6 cm h" 1 ; for this last value,' the input power 
has decreased to 14.05 kW. From Fig. 20, we find that 
the maximum pulling rate is 4.4 cm h" 1 since, at 4.6 
cm h" ! , the melting point isotherm extends within the 
melt around the crystal ; it is therefore impossible to 
. rn a in tain a crystal of constant radius. Even at 4.4 cm 



Table 1 . Thermal properties of the materials used in the germanium furnace 





Symbol in 


Conductivity 




Material 


Fig. 18 


(Wm-'K- 1 ) 


Emissivity 


Solid germanium 


GS 


25.0 


0.55 


Liquid germanium 


GL 


75.0 


0.20 


Graphite 


G 


60.0 


0.81 


Steel 


S 


30.0 


0.30 


Insulating material 


F .. 


0.07 
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Fig. 19. Sketch of the germanium furnace with finite element 
mesh (left) and isotherms separated by steps of 10 K (right). 



tr \ we find that the melting point isotherm is tangent 
to the surface of the melt; such a temperature dis- 
tribution could make it difficult to control the diam- 
eter of the crystal. Quite clearly , the effect of increasing 
the pulling rate is to produce a concavity of the crystal 
at the interface, with a deflection of 8:6 mm at a 
growth rate of 4.4 cm h" 1 . To the contrary, a low 
pulling rate of 1.4 cm h" 1 leads to a convex crystal,' 
with a deflection of 6.9 mm, as we may see in Fig 21 • 
the input power is then 25.0 kW/ '" ■ ? 




Fig. il. Isotherms in the' germanium jfurnace separated by 
steps of 10 K for an imposed growth rate of 1.4 cm h" 1 . 



It is clear that slight- changes of the input power 
(between 24 and 25 kW) are associated with major 
modifications of the. pulling rate (between 4.6 and 1 A 
cm h ') for keeping the crystal diameter constant ; it 
was therefore' advisable to use the pulling rate rather 
than/the power input as the-control parameter for the 
calculations. Let us however analyse the case -where 
the power input is imposed at a value of 26 kW, with 
the result at 25 kW as a first guess. With such a jump 
in the input power, it is impossible to : apply in a 
straightforward manner the method of Section '6, 
because the initial guess lies too far from the solution. 
However, convergence can be obtained by applying a 
sub-relaxatioh algorithm, i.e. by imposing as con^ 
straint on the interface. Instead of equation (80), we 
use ; • 

7^ +1 .=.(i-a)r m +ar U ' (88) 
where ? is a sub-relaxation, parameter while n refers 
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Fig. 20. Isotherms separated by steps of 10 K in the germanium melt and crystal, for'different growth 
rates:(a)4.2cmh- 1 ;(b)4.4cmh" , ;(c)4.6cmh- 1 . 
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Fig. 22. Evolution .of growth rate (a), temperature at the tri-junction (b) and position of the interface 
along the axis of symmetry (c) during the numerical iterations when the power input is imposed. 



to the nth global iteration. We found that an initial 
value of 0.75 for a was suitable, in this case, for con- 
trolling the oscillatory behaviour of the iterations. 
Figure 22 shows the evolution of the pulling rate t? pnl 
(a), the imposed interface temperature T x (b) (see 
equation (78)), and the axial interface position (c) dur- 
ing successive iterations. Having started the procedure 
with a = 0.75, we select the value 0 which irnmediately 
generates oscillations. A strong damping is observed 
once the value of 0.5 is selected. The final solution is 
shown in Fig. 23 : the crystal becomes much more 
convex at the interface (with a deflection of 22 mm) 
while the pulling rate reaches a negative value of —2.2 
cm h~ l . Since the power input is too high, we are actu- 
ally melting the crystal by pushing it into the melt : the 
natural issue would be to grow a crystal of smaller 
radius. 

8.2. Gallium arsenide growth 

As a second example, we will calculate the LEC 
growth of a gallium arsenide crystal. The liquid boric 
oxide layer, which is used to prevent evaporation of 
the volatile component, can be considered as trans- 
parent for a range of wavelengths and opaque for the 
others. It has been shown in ref. [43] that transparency 
occurs for wavelengths lower than about 2.3 fim. 

As explained in Section 3.4, we will consider two 
different radiative enclosures, corresponding to both 
ranges of wavelengths. For the first enclosure, the 
function y A (T) defined by equations (2), (3) and (4) is 
given by 



y A| =F(2.3* 1511) -f(0) = 0.4 (89) 
while for the second enclosure, we have 

7a 3 = F(oo) -F(2.3 * 1511) = 0.6. (90) 

We wish to investigate the effect of the value of y Al 
upon the global solution. Recall that y Aj = 0 for a 
fully transparent medium, while y Aj = 1 when it is 
fully opaque. «~ : 

The physical properties adopted for the present 
simulation are given in Table 2. Note that the melting 
temperature of gallium arsenide is 1511 K. The 
geometry of the furnace and the finite element mesh 
used for the simulation are shown in Fig. 1 . The crystal 
diameter is 6.4 cm (2.5 in.), the length is 5.2 cm and 
the total weight (melt and crystal) is 4 kg'. 

A global view of the isotherms for the case y Ai = 0.5 
and a pulling rate of 1 cm h~ 1 is given in Fig. 24. We 
note that, for the present problem, we have calculated 
the shape of the meniscus for the surface of the melt 
and of the oxide layer. 

In Fig. 25, we show the isotherms and the contour 
lines of the Mises invariant 5 M defined by equation 
(87), for values of y^ respectively equal to 0, 0.25, 
0.50, 0.75 and 1, and for the. same pulling rate of 1 cm 
h~ l . It is clear that the degree of opacity of the oxide 
layer has a major impact upon the process. In particu- 
lar, in going from the transparent case at y A2 = 0 up 
. to the opaque situation at y Ai = 1, we find that : 

(i) the crystal, which is initially concave at the inter- 
face, becomes convex ; 



Table 2. Thermal properties of the materials used in the gallium arsenide 

furnace 







Conductivity . 




Material 


Symbol 


(Win" 1 K- 1 ) 


Emissivity 


Solid gallium arsenide 


Ga As S 


1J2 


0.36 


Liquid gallium arsenide 


Ga As L 


17.1 


0.36 


Graphite 


G 


42.0 


0.64 


Quartz 


Q 


3.0 


0.5 


Steel 


S 


• 27.2 


0.45 


Boric oxide 


B 


1.7 


0.5 


Graphite felt 


F 


1.0 
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Fig.' 23. Isotherms in the -germanium furnace separated by 1 
steps of 10 K for an imposed power input of 26 kW/ . 



(ii) the vertical temperature gradient in the melt 
decreases, while the horizontal temperature gradient 
on the surface of the melt increases ; * 

(iii) the vertical temperature gradient in the oxide 
layer strongly increases ; 



1300- 




Fig. 24. Isotherms in the LEG furnace for yx - « 0.5 and 
pulling rate of 1 cm h" 1 , separated by steps of 50.K. 
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(iv) the stress level increases along the edge of the 
crystal ; 

(y) the stress level along the axis of symmetry 
decreases near the interface. 

The sensitivity of . the system is summarized in Fig. 
26, showing the dependence of some important output 
parameters upon y A2 , i.e. the maximum value of S M , 
the power input, the deflection of the interface and 
the total heat flux leaving the upper surface of the 
boric oxide layer. We note that all these parameters 
depend appreciably upon the radiative properties of 
the encapsulant. 

9. CONCLUSIONS 

The accurate prediction of the shape of the melt- 
crystal interface and the temperature field within the 
crystal requires a global modelling of heat transfer in 
the entire furnace. We have presented a self-contained 
model where, for any given geometry, the input par- 
ameters are limited to the coolant temperature, the 
power input or the pulling rate, and the radius of the 
crystal. Radiative exchanges are accurately calculated 
since viewed and hidden parts of the enclosures are 
taken into account. The effect of melt convection will 
be analysed in a further publication. 

We have demonstrated, on the basis of simulations 
of germanium and gallium arsenide growth, that our 
model is relevant for evaluating the sensitivity of the 



final product upon material and geometrical par- 
ameters. It is a valuable tool for growth control and 
furnace design. 
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Global modelling of heat transfer in crystal growth furnaces 

UN MODELE GLOBAL POUR CALCULER LES TRANSFERTS DE CHALEUR DANS 
LES FOURS DE TIRAGE DE CRISTAUX v ^ 

Resume — Pour predire en cours de croissance le champ de temperature dans un cristal, ainsi que la forme 
de r interface qui separe le bain de ce cristal, il faut determiner avec precision les transfsrts de chaleur qui 
ont lieu dans le four tout entier. La resolution de ce probleme est tres complexe, car il y a lieu de calculer 
precisement a la tois les echanges radiatifs entre les difierentes surfaces et les transferts par conduction a 
rinterieur des constituants du four. -Les echanges radiatifs sont evalues sur base de l'hypothese que la 
radiation est diffuse et se fait en surface tandis que les calculs se font par la methode de discretisation de 
Galerkin, a Taide d*un algorithme specifique de calcul des vues et cachees. Le modele a ete etendu pour 
prendre en compte le caractere semi-transparent de certains materiaux. La forme de l'interface liquide- 
solide est une inconnue du probleme, que Ton calcule en Fassimilant a Tisotherme de fusion. Des exemples 
de fours de tirage de germanium et d'arseniure de gallium sont analyses, ce qui permet d'illustrer la 

puissance de la methode. 



EIN GLOBALES MODELL FttR DEN WARMETRANSPORT IN OFEN FOR DIE 

KRISTALLZOCHTUNG 

Zusaramenfassung — Eine quantitative Berechnuhg des Temperaturfeldes und der Position der Kristall- 
oberfiache beim Wachstum erfordert eine genaue Kenntnis der Warmetransportvorgange im gesamten Ofen. 
Dieses Problem ist auBerst komplex, da es eine exakte Berechnung der Strahlung zwischen den yer- 
schiedenen Oberflachen und der Wanneleitung in den einzelnen Komponenten erfordert. Der Strahl- 
ungswarmeaustausch wird unter Annahme von diffusen grauen Oberflachen berechnet. Hierzu wird ein 
besonderer Algorithmus zusammen mit einer Diskretisierung nach Galerkin verwendet. Dieses Modell 
wird erweitert, urn auch halbdurchlassige Materialien einbeziehen zu konnen. Die Form der Flussigkeits/ 
FeststofF-Grenzflache ist eine Problemvariable — sie wird in Form der Schmelzisothermen berechnet Die 
Wirksamkeit des Verfahrens wird anhand von Beispielen fur Germanium- und Gallium-Arsenid-Ofen 

gezeigt. 



MOAEJIHPOBAHHE TEIUIOIIEPEHOCA ITPH POCTE KPHCTAJDIOB B EMKOCTH 

AHHOTaqHs — Rjis. KOJiHHecTBeHHoro onpeaejieHHH TennoBoro nojia h noJioxeHHH rpammH pa3aejia 
pacruiaB-KpHCTajui b nponecce pocra Heo6xonHMu TOHHue gamme no TemionepeHocy bo bcch 
eMKocTE. 3ra 3aaana BectMa cjioaoia, tuk Kax cB3ana c tohhhm pacneroM nepeHoca B3jryHeaB5i Mexyry 
paanEHHLiMH noBepxHocTSMH h TenjionpOBOAHocTbxo bo bccx hrctzx CHcrreMBL JiyHCTBrfi tchjiooomch 
paccHHTbmaeTca b ripHOJiHxeHHH m*4><|>y3Ho-cepwx noBepxHOXTefi h c HcnojiB30BaHHeM ajrropHTMa 
BHflHMofi e CKpbiTOH ^acTCH, a TaicKe Meroaa TajxepcHHa. Modern, o6o6ineHa Ha yner nojrynpo3pa*iHLix 
MaTepnajioB. <X>opMa rpammu pa3,nejia saaztKocrb— TBep^oe tcjio sBJiaeTCS nepeMeHHofi BejnranHOH 
3TOH 3ajnaHH h onpcaejiaercji KaK rooTepMa ruiaBjieHHs. Asajna npoueccoB rjik apcemmoa repMamui h 
rajuiHH inoffTBepOTJi 34>4«kthbhocti» MeTona. 
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